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TBCWRARY STORAGE 

PARAMETERS 
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DURATION 



DESCRIPTION 



Solution of a System of Differential Equations "by Milne's 

Iterative Method (SADOI Only) 

Closed subroutine 

Interlude 80 Subroutine h6 

Interlude 0, 1, 2 Subroutine 0, 1 

13 locations specified "by S6 

S^ - 88 

These parameters «*"d the computed parameter stored in 

location 9 should not he destroyed during the subroutine. 



00 F 00 aF 



a is location of first word of auxiliary 
routine 



00 F 00 hF N(h+i) are the independent variahles x 



00 F 00 cF 



7 00 F 00 nF 



8 00 F 00 hJ 



Time of Interludes I(ll«3 +3 2 

i=0 

Trp^imm (i ) + .5 milliseconds 



oi 
(i = 0, 1,.2, .<., n-l) 

Locations c thru c+12 are used as 

temporary storage 

n is the number of differential equations 

to he solved 

h is the length of each step of integration 

n-l 



D 1 ) n - .1] x 



Time of Subroutine: 



n-l 
I 2.5 I + (I.+1) D -i- 13 -5n 

i=0 1 1 



milliseconds per step of integration 

n - number of equations 

I = number of iterations required per equation 

D = duration of auxiliary routine 
This routine will solve n simultaneous first order ordinary 
differential equations expressed explicitly as a function 
of the variia&Les. It is also possible to solve an equation 
of order n "by reducing the equation to a set of n first 
-order differential equations. 



-a* 



EXAMPLE 



PROCEDURE 



y" + I y' + o g(x) » can "be written as 



*0 



y i * "* y i " * s ^ 



At location h, specified "by a parameter in location 5, 

oi 
y" , y 1 , and y . Before the program can be hegun it 



oi 



oi 



oi 



is necessary to have three other values of the function 
and their derivatives « These are all f ound by an interlude 
using several equations devised hy V. E. Milne* The infor- 
mation is stored at the following locations? 



~S 



I (d + 1 + % ) = y_ 



11 



(i « l, 2, *.*, n-1) 



II (ft + i + 



- v-t 



li 



H (b + 1 + 5n) = y. 



11 



I (k + I + 6n) = y 



21 
I (b + i + Tft) * y 21 

After these values have heen found, the subroutine Is 
automatically read in, and the main routine can fee begun* 
Each time the subroutine is called into use it will carry 
out erne integration step of length h and store the value 
of the function. The subroutine will also find the value 
of the fieri vat ive and store it. These values will "be stored 
in the following manner! 

y* « BTb + 2(k + l) n + 1] 

k = 0, 1, . .., M where ,M is the number of 

integration stejs 

i = 0, 1, ••«, n-1 

y^ = Ifc + (2k + 5) n + i] 

k = — 1/ 0, lj 2, ***jr M 
i * 0, 1, 2, *.«, n-1 
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During the interlude the value of x oscillates from x Q + h 

to x - h and then changes to X- + 2h. Therefore it is 

^ 

necessary for the coder to know -where x is at all times . 

For each equation x can he found at 1S6. 

Milne's method uses two quadrature formulas which 

first t&*edict a value of the function and then make a 

correction. When the correction is less than or equal 

to 2"^, then the new value will he accepted. Since the 

ptfbgrammer does not know exactly how many iterations are 

necessary each time, a counter hfcs "been provided which will 

tell the coder which equation is "being used. In location 

-39 
S6 5 during "both the interlude and subroutine i x *2 

•will he found. 

1$ is jpssihle that h has "been chosen too large so 
that the predicted value will he of little use. If this 
Is tru$, then a carriage return and line feed character 
will he punched followed hy several F*s m& then the machine 
will he irtopped. Dividing h hjr 2 will decrease the error 
%j a factbr of* thirty-two. 

During the entire operation the ahsolute value of 
%h@ derivatiTe must he less than one-half. If the derivatives 
are greatdi than or equal to one-half then overflow will 
occur and this will also cause F 1 * to he punched. Finally, 
the punching of F's may also indicate that 4 faulty sub- 
routine has "been written. Before using this code, the 
programmer should determine hy usi^g a code check the average 
number of iterations necessary. Of course, the number of 
iterations wfl.ll vary with the diff lenity of the equation. 
DlSCUBSIOfl" far MftTHEMSTICAL METHOD USED 

A. Int erlude 

m*mm •«»*.»«— ^ hot m 

Trial values of y' and y^ are computed from the relations 



JW 



These -values are substituted into equations JL and iL 

to giye first approximations to y- and y , • The first 

approximations are then substituted into the given differential 

equation and improved values of y' and y* are computed. 

The process is continued until the change in y' and y' 

-39 - 1 ** L 

is not more than + 2 • 

The next step is to use equations B and C in a similar 

-39 
way to obtain a value of y 2 , the tolerance again being + 2 « 

Here ve ohtain y from B, y' from the differential equation 

ftT>d j again from C» 

(kj y x = y Q + h/2* (y^ + 16 y * Q + 7 y{) + yj h 2 A 

c^) y-i - y fc h/2Jf (7 y il + l6 y + *V + ^ h2 / 4 

(B) y 2 - y Q + 2h/3 (5.y£ - ar£ - y^) - 2 yj h 2 

(C) y 2 = y Q + n/3 CyJ + ^ y£ + y£) (Simpson's Rule) 

B. Subroutine 

With the initial values and the three values found toy 

the interlude an approximation y^ * to y.. is made toy using 

formula (B). This value £s then substituted into the 

differential equation and the result is use£ in formula 

(2) 
(E) to obtain an approximation y^ 7 to y,. It should 

fee noted that (s) is simply* Simpson's rule* JUT y^ ' 

and y^ ' agree, then we are finished} otherwise, vith 

y.jj 2 ' an Improved evaluation of y* is made and the process 

is continued. 

«»> y »*l = y »-3 + te/5 (2y i-fi " y i-l + ^ 

< E > y *n = y *-i + ^3 (y;.i + V; + y^ x ) 

It can he shown that the final correction id a factor 
of the original correction C Q . That is, C m * finC Q vhere 
© = hy'/3. 
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Tlie process will converge if and only if |©| < 1. 
Dividing h "by 2 will not necessarily reduce the iterations 
by a factor of two. The number of iterations per step 
of integrations will "be reduced "by a factor of two or more 
only if > 1/2 . 

It is now possible to estimate the approximate number 

of iterations necessary per step of integration- The final 

-39 
correction is never greater than 2 and the original 

-12 
correction is never greater than 2 or else the machine 

will stop and punch an P. 
C. Error 

The truncation error can be estimated by the computation* 
_ „(!) j. is- . -n „ Pfl/on y?Tr&)~ 



Equation (D) y = y££ + ^ , \ = 28/90 h P y 

r (2) 

n+1 



(E) 7 = ?{ 2 J '+ E c * E ■ = - h 5 /90 y (5) 



y n+l + y n+l = " 29h5 / 90 y (5) = 29 ( * ^ 9 ° h5 y(5)) = 29 \ 

*2 v 'n+l J n+1' ' 

It is also possible to make an estimation of the maximum 
accumulated error of a numerical integration over a range 
of length L with N ec|ual steps. Let G be a positive constant 
such that |y n | < G an£ let M be a positive constant such 
that \y^) | < M , then E^ < l\ Je 210 " 1 )/l80 N* G. 

The value of G t ftn*t be estimated from the computational 

That is, G - (y£ - 7^ m2 )/(y n " y n ^) * ? can " be art*"** 8 * 
by knowledge of the trueation error* Since Eg = (-h 5 y^')/90» 
y ^) = (^o E )/h . Uaturally, the error of the subroutine 
will be dependent on the error of the values found by the 
interlude. For complete analysis of the error, see Richter, V. 
"Sur l'Erreur commisse dans la methode d' integration de Milne", 
Comptes Rendus de l'Acp&emie des Sciences, vol. 233 (1951) 
pp. 13^2 - 13^. 
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LOCATION 



CRDER 



00 K(J2) 
26 1000H 






00 F 




00 3F 


1 


1*5 1S6 




40 S5 


2 


00 F 




00 L 


3 


38 P 




00 F 


k 


51 8F 




66 L 


5 


S5 S7 




kO 9F 


6 


4l S6 




4l F 


7 


L5 S5 




40 1S6 


8 


19 38F 




lA- 7L 


9 


42 7L 




IA 5L 


10 


46 7L 




19 IF 


11 


Ik F 




40 F 


32 


36 7L 




50 8f 


13 


75 2S6 




40 F 


Ik 


IA 5S6 




40 5S6 


15 


LI F 




IA 3S6 


16 


40 10S6 




50 F 



NOTES 
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F 2 



7/16 



from 76L 



h/3-»9F 



Bring out * 01 , y^, y^, and y Q± and 
store at 1S6, 2S6, 3S6* 4s6. 



Count 



y o h 



rf o ^o 



rf o ^o 



LOCATIOH 


CEDES 


17 


V 8f 




10 2F 


18 


^0 966 




L5 1S6 


19 


IA 8F 




1*0 1S6 


20 


k-S ns^ 




L5 10S6 


21 


10 Iff 




IA 3S6 


22 


40 F 




50 5S6 


23 


7J 3L 




lA F 


2* 


ho 7 




50 F 


25 


75 9F 




00 IF 


26 


Ik 9S6 




IA 4s6 


27 


40 (6s6) 




50 27L 


28 


26 Sk 




kQ F 


29 


LO 5S6 




kO IF 


50 


19 38F 




12 IF 


31 


36 32L 




4l 11S6 


32 


L5 10S6 


i 


40 5S6 


33 


L5 F 




40 10S6 



by k2, kk 



I&fES 



E&GE 2 



y 



h 2 A 



x + h 

o 



Store l/2 at 11S6 

V ^f &* + 16 y; + t y£) 



Form 



(y n h*" 



)A 



Call in auxiliary s*dbrpufc5:&e 



Test 2 
if < 



39 



- \r™-riW 



r 



n 



J 1 



Ki 



j**j±. 



r x aBd r x are laterchanged alternatively. 



LOCATION 


ORDER 


NOTES 
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34 


4l 12S6 
50 12S6 








35 


LI 8F 
40 8F 








36 


00 IF 
L4 1S6 








37 


40 1S6 
LI 9F 




x + h 

— 




38 


40 9F 
LI L 




± V3 




5? 


40 L 
32 41L 




Binary switch 




40 


L5 13L 
46 27L 








in 


22 20L 
L5 79L 








42 


46 27L 
L3 11S6 




Test for zero 




43 


36 20L 
L5 1S6 








44 


L4 8F 
4o 1S6 




X' + 2h 






45 


L5 3S6 
L4 10S6 








46 


10 3F 
kO F 








47 


L5 5S6 
10 IF 








48 


40 IF 




Fora y Q + (2h)/3 (5 7[ - 7* Q - 7\±) '■- ^7 n Q ** 




10 2F 








*9 


L4 IF 
LO F 








50 


40 F 
50 F 



















LOCATION 



Pi 



52 
53 
54 
55 
56 
St 

59 
60 

fit 

& 

65 

66 



j 



ORDER 

75 9F 
00 IF 
LO 9S6 
00 3F 
l4 4s6 
4o 8s6 
22 5^L 
50 5^L 
26 S4 
40 ts6 
L4 3S6 
10 2F 
L4 5S6 
40 F 
50 F 

75 9F 

00 2F 
L4 4s6 
40 F 
LO 8S6 
40 IF 
L5 F 
40 8s6 
19 38F 
LO IF 
36 65L 
L5 8S6 
22 54L 
L5 1S6 
40 S5 
19 18F 
L4 65L 
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from 64 



Auxiliary routine 



Form 



y + */3 (y; + ^{^ y£) 



(2) „ (1) 



- X 



a" 39 - | y | 2) - y, 

if > 
if < 



(i) 



Store values permanently 



LOCATION 
6? 

68 

69 

70 

71 
72 
73 

7* 

75 
76 
77 
78 
79 



ORDER 

L4 IF 
ho 6% 

19 2F 
IA 12S6 
UO 12S6 
36 65L 
19 38F 
l4 S6 
40 s6 
LO 7F 

32 76L 

L5 S6 
L4 1L 
40 65L 
00 2QF 
46 7L 
L5 36L 
42 7L 
22 6L 

L5 5F 
L4 IF 

40 5P 
50 2L 
26 999F 
00 6s6 

26 4L 
26 1H 

S5 F 
L4 l4L 
42 44L 

41 s6 

L5 S5 
40 1S6 

19 5F 
4o 986 



HOTES 



1AGE 5 



F 2 



Does 1 = at 

If I -i n change addresses in 65 and 66 



from 72 



Find location of y, and store at 3 



Set link address 



Bring out xl. 



LOCUTION 
4 

5 

6 

T 
8 

9 
10 
11 
12 

13 
14 

15 
16 

17 
18 

19 
20 
21 



ORDER 

L5 S3 
40 2S6 
19 38F 
L4 4l 
42 4L 
L4 32L 
46 4L 
L5 9S6 
L4 9S6 
4o 9S6 
36 4l 
L5 4L 
10 2QF 
l4 7* 
^2 33L 
L4 7F 
42 31L 
LI 5S6 
10 IF 
L4 3S6 
L4 7S6 
4o IF 
50 IF 
75 9F 
00 3F 
L4 2S6 
40 8s6 
50 i7E'.- } 
26 S4 
L4 5S6 
10 2F 
L4 786 
40 F 
50 F 

75 9* 
00 2F 



NOTES 



KAGE 6 



F 2 



Bring out initial values and store 



Count 



Set store address for y* . 

ki 



and y 



ki 



Form 



y n . 3 + (V5) (2 jr^, - y ;. x + 2 y ;) 



Call In auxiliary subroutine 



Form y - + (h/3) (y' ,+^y' +y' J 



locatioi 


ORDER 


IOTHB PAGE 7 




22 


IA 6s6 


| 






40 IF 








23 
2k 


LO 8S6 
kO F 
L5 9S6 
32 28L 




Store y^ - jl 1 ) 
J n+1 'n+1 

Test for first iteration 




25 


19 l^F 
12 j 








26 


36 28% 
92 I29F 




is |A> . y (l) |>2- 12 
Ho, transfer to 28 and continue 




27 


92 898F 
OFF 




Yes, print F's and stop 




28 
29 


41 9S6 
19 58F 
12 F 
36 31L 




Is | v ( 2 ) . yWl < P -39 
Yes, transfer to 31 




30 


L5 IF 
26 17L 




Ho, transfer to 17 




31 


L5 IF 










ho ( )F 


>y 12 


Store y n+1 




32 


50 S7 
50 32L 




' a 




33 


26 Sk 




Call in auxiliary subrctttine 




3* 


ho ( )F 
19 38F 
IA 86 


T?y U 


Store y 1 „ 




35 


40 S6 
LO 7F 




Increase i 

i - n? 




36 


32 kOL 
51 S6 




Yes, transfer to kO 

Ho, increase addresses by 1 and continue 




37 


00 59F 
Lfc 4% 




to solve other differential equations 




38 


40 4l 
19 18F 








39 


IA 2L 
40 2L 




... 





7 2 



LQCAHnOH 

*5 



ORDER 
26 2L 

00 60F 
I* **5L 
1*0 kl 

L5 khl 
kS 2L 
50 85 
22 ( )F 
L5 S3 
1*0 2S6 



HOTES 



BWS 8 



from 56 



Increase addresses by 2n 



by 1 



Constant 

Leave subroubiae 
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